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In this paper we introduce a continuous time stochastic neurite branching model closely 
related to the discrete time stochastic BES-model. The discrete time BES-model is underly- 
ing current attempts to simulate cortical development, but is difficult to analyze. The new 
continuous time formulation facilitates analytical treatment thus allowing us to examine the 
structure of the model more closely. We derive explicit expressions for the time dependent 
probabilities p(7, t) for finding a tree 7 at time t, valid for arbitrary continuous time branch- 
ing models with tree and segment dependent branching rates. We show, for the specific 
case of the continuous time BES-model, that as expected from our model formulation, the 
,-0 , sums needed to evaluate expectation values of functions of the terminal segment number 

fi(f(n),t) do not depend on the distribution of the total branching probability over the ter- 
minal segments. In addition, we derive a system of differential equations for the probabilities 
p(n, t) of finding n terminal segments at time t. For the continuous BES-model, this sys- 
tem of differential equations gives direct numerical access to functions only depending on 
the number of terminal segments, and we use this to evaluate the development of the mean 
and standard deviation of the number of terminal segments at a time t. For comparison we 
discuss two cases where mean and variance of the number of terminal segments are exactly 
solvable. Then we discuss the numerical evaluation of the S-dependence of the solutions for 
the continuous time BES-model. The numerical results show clearly that higher S values, 
' i.e. values such that more proximal terminal segments have higher branching rates than 

more distal terminal segments, lead to more symmetrical trees as measured by three tree 
symmetry indicators. 
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I. INTRODUCTION 

Trees appear in many areas of science. If we limit ourselves to some examples of naturally 
occurring tree-like structures we find streams (rivers, creeks) in geology (e.g. Horton [2l|], Shreve 
[3c]]). actual trees and shrubs in botany (e.g. Bell et al. {5], de Reffye and Houllier 9j, Godin 



et al. 
et al. 



14], Sismilich et al. 



111 ]. Dityatev et al. 



3111 ) and axon and dendrites in neuroanatomy (e.g. Asco 
12] ]. van Pelt and Schierwagen [^j], Verwer and van Pelt 



i I If] , Devaud 



451] ). For our 



purpose, which includes the description of both the temporal development of the tree population 



3 



and the resulting final population, there is an important distinction between streams on the one 
hand and plants and neurites on the other hand. River networks are the result of a merging of 
streams originating from independent sources. Plants and neurites, however, grow from the root 
and are the result of branching and pruning. We develop our formalism for the latter, i.e. processes 
where tree-like structures grow from an initial root-like element. 

Modeling dendritic and axonal morphology is a field of growing interest in computational neu- 
roscience and has led to the formulation of several types of models: statistical population models 



(Ascoli Nowakowski et al. [27], Samsonovich and Ascoli 29|, van Veen and van Pelt 43]), Lin 



denmayer or L-systems (Ascoli and Krichmar [2|, Torben-Nielsen et al. 32]), biophysical models 



(Graham and van Ooyen [15J, Hely et al. [18], Hentschel and Fin e |19l| . Hentschel and Van Ooyen 
201 ]. Kiddie et al. [23J), and stochastic growth models (Kliemann [251 ] . Uemura et al. 33J, van Pelt 



et al. 



41]). These models and their extensions are also underlying current attempts at large scale 



modeling of cortical structures (Eberhard et al. [13], Koene et al. 



26(], Zubler and Douglas 48]). 



Exact solutions provide a solid reference against which aspects of the algorithms generating the 
neuronal morphologies for these large scale network models can be tested. 

The BES-model, whose name derives from the symbols used for its main parameters, was 
proposed by van Pelt et al. 



41 



421 ] and is the main topic of this paper. The BES-model was 
developed as a stochastic description of the process of neurite outgrowth. The goal of the BES- 
model is to capture not only the composition of the final population of trees but also its temporal 
development. Here we develop the mathematical and computational tools for the evaluation of this 
and more general models of branching trees and try to elucidate the structure of these models. The 
BES-model is a conceptual merger between the BE-model describing the temporal development 
of the number of terminal segments and the S-model describing the competition between between 
terminal segments. It has been shown by Villacorta et al. [46] that the outcomes of the BES-model 
cannot be factorized into an S-model and a BE-model contribution, thereby challenging the idea 
that the BES-model is a valid combination of the two. The reason this factorization fails probably 
lies in the fact that multiple branching events can take place in a single time step. In our treatment 
of the BES-model, in which we reformulate the model in continuous time, we will show that the 
dependency on the B and E parameter can be separated from the dependency on the S parameter 
and that therefore a factorization into a BE-model and an S-model is actually possible. This leads 
us to the very generic idea of ir, p- models in which ir stands for the dependence of branching rates 
on the location of the terminal segment in the tree and where p captures the overall branching 
probability of a tree. 
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In biology branching processes are often studied for the description of population dynamics. 
This is a tradition which dates back to the early work of Galton and Watson on the survival of 
family names. There is a natural link between branching processes and trees as the description 
of the history of a branching process constitutes a tree. The theory describing the composition 
and number of terminal segments in branching processes is well developed (Asmussen and Hering 



,3, 



, Athreya and Ney [4], Kimmel and Axelrod 



24]). Unfortunately for our purposes this focus on 



the composition and number of terminal segments disregards the history of the individual branching 
process. As a result relations between parent and children are only used to describe the passage 
from one generation to the other. For branching biological processes like dendrites or axons and 
also for real trees the current shape of the tree is a direct reflection of its history and the different 
branching events it experienced. Hence keeping track of the individual histories becomes important 
and the main objective of this paper is to show how this can be achieved correctly and efficiently in 
those cases where only the terminal segments of the tree can branch. The more general n, p-model 
which we use in our analysis can include intermediate segment branching as well. 

We develop the theory for the continuous time BES-model from that of the Poisson process 
(Cox and Lewis Dehling and Kalma [3]) by including branching into the Poisson process. The 
terminology 'branching Poisson process' has been used before to describe interactions between two 
or more Poisson processes in which events in one process influence the Poisson rate in the other 
process. These Hawkes branching point processes (Hawkes are suited to model the interaction 
jetween firing rates in different populations of cells (Cardanobile and Rotter Chornoboy et al. 
3] , Johnson [22] ) or the influence of a large earth quake on the incidence rates of small earth quakes 
and vice versa (Veen and Schoenberg 44], Zhuang 43]). Hawkes branching point processes deal 
with a limited set of mutually interacting Poisson processes. In tree growth, on the other hand, 
every branching event is a transition from one stochastic process to another and not a repetition in 
a chain of equivalent events. The important shared aspect between the theory of Poisson processes 
and the theory of branching trees in continuous time are the probability density functions used to 
describe branching and survival of tree segments. 
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II. MODEL AND ANALYSIS 
A. The BES-model in discrete and continuous time and the tt, p-model 

In the topological view on branching which we use here a segment is defined purely on the 
basis of the topological properties, i.e. a segment is part of the tree that connects two branch 
points (intermediate segments), or it connects the root point to the first branch point (initial or 
root segment) or a branch point to a terminal point (terminal segment). We limit our treatment 
to binary trees, for these trees given a total number of terminal segments n the total number of 
segments including intermediate and root segments is 2n — 1. When we evaluate the mean of 
the number of terminal segments and and its variance, this choice of n limits our need for extra 
multiplication factors and additions and thereby simplifies our calculations. 



As pointed out before the original BES-model was proposed van Pelt et al. 
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42] stochas- 



tic model for describing changes in dendritic topology during neurite outgrowth and has been 



successfully applied to categorizing morphological data (Dityatev et al. 12j], van Pelt et al. 42j]). 
The BES-model uses the parameters B,E, and S to set the branch probability p sn per time step 
At of an individual terminal segment s in a binary tree 7. Let us introduce these parameters. The 
parameter B sets the probability of branching for a tree with only one single terminal segment. 
The other two parameters are used to relate the branching probability of a terminal segment s 
in an arbitrary tree 7 to B while intermediate segments don't show branching in the BES-model. 
The overall probability p 7 of a tree 7 to branch depends on the number of terminal segments ra 7 
and is given by 



Pi 



Bn^ E \ (1) 



This probability is distributed over the different terminal segments s depending on their centrifugal 
order v s , where the centrifugal order counts the number of segments separating the tip of a segment 
from the root point. Given that a tree 7 branches, the probability ir s that this will take place at a 
terminal segment s is given by: 

2~ Su a 



with C y = J2 2 ~ Sua - ( 2 ) 



^7 SG7 



Combining tt and p we obtain for the branch probability p Sj7 : 

o—Sv, 
1 ^ 



p S)1 = p^TTs = Brtf E) — — , (3) 



6 



with, because we are dealing with probabilities, /o 7 ,7r s > for all trees and segments. We have to 
point out here that we, contrary to the original formulation of van Pelt, have moved the terminal 
number dependence out of the normalization factor and made it explicit as an extra contribution 
to the exponent of re, i.e. we write n 1 ~ E with C 7 , where van Pelt writes n~ E , and uses C 7 = ^-C 7 
as a normalization constant. Our choice turns the centrifugal order dependence into a simple 
distribution of the total branching probability over the terminal segments, i.e. X)sG7 7r s = 1- As 
will become clear in this paper, our choice better captures the real structure of the problem. 

We will call the BES-model as described above the discrete time BES-model (dBES-model), 
because it can be seen as discretized version of the continuous time BES-model (cBES-model) which 
we will introduce shortly. Our original transfer operator-based analysis of the dBES-model (see 
appendix[B]) made us aware of the limitations of the discrete time formulation. In the discrete time 
formulation the counting of different branching histories leading to a specific tree is complicated 
by the presence of time steps without branching events. In addition the discrete time formulation 
introduces a small S'-dependence in the total branching probability of a single tree because it 
affects the probability of finding more than one branching event in a time step. We then realized 
that counting is much simpler in continuous time because we can limit ourselves to counting the 
different orders in which terminal segments branch while the exact timing can be dealt with by 
taking integrals over probability densities. It was this realization that counting histories in the 
continuous time would be much simpler that led us to formulate and explore the continuous time 
BES-model. 

In the cBES-model tree growth is modeled as a branching Poisson process. Poisson processes are 
characterized by rates instead of probabilities, and consequently the continuous time BES-model 
specifies a branching rate for every terminal segment of a tree. This branching rate applies until 
the next branching event takes place after which the new terminal segments and the surviving 
terminal segments get new branching rates. The branch rates A Sj7 are chosen in such a way that 
to first order in the time step we get agreement with the discrete model, 

X s = 2m. = bn\~ E ^— with b = ±L. (4) 

' 7 At 7 C 7 At w 

The probability p(^y, s, At) that only the Poisson process associated with the terminal segment s 
in the tree 7 produced a single branching event can be obtained by integrating the corresponding 
probability density. This probability density is a product of the probability densities for the 
branching A Si7 e _As ' 7 * of segment s and the survival fls'^s e _As ' 7 * of segments other then s, and we 
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find that: 



p(j,s,At) = / A S)1 e- A "'n^^4 

->0 77„ 



s'/s 
/■Ai 

= ir s p 1 e~ p ' yt dt 
Jo 

= tt s (1 - e^ At ) 

= 7r s (p 7 At + h.o(At)) 

= p s + h.o(At). (5) 

There is a subtlety in this mapping which might easily go unnoticed. The dBES-model allows 
synchronous branching of multiple segments, but with probabilities which are higher order in At. 
These higher order contributions are kept low by keeping B or equivalently At small. We, however, 
calculated p Sjl from A Si7 under the assumption that only the segment s branches. The cBES-model 
allows for an analysis in which only direct transitions between trees differing by one branch event 
take place, while analysis of the dBES-model needs to incorporate transitions between trees differing 
by more than one branch event. 

In a large part of our treatment of the cBES-model we only use part of the structure obtained. 
To keep our results as general as possible it will therefore be advantageous to choose a notation 
which fits this structure and which leads to a more abstract model, not limited to the cBES-model. 
In this tt s , /> 7 -model (or 7r, p- model for short) the whole tree branching rate p 1 depends solely on 
the tree and is then distributed over all terminal segments, 

As, 7 = TtsP^j 

J>. = 1- (6) 

Above we used the notation s £ 7 to indicate that s is a terminal segment in the tree 7. Furthermore 
because each A Sj7 represents a rate we take each p 7 > and each ir s > 0. The cBES-model is a 
subclass of this set of models as is clear from the following identifications: 

2— Sv s 



p 7 = bn 1 . (7) 

We limited the 7r, p-model to terminal branching, but most results derived here apply equally well 
to non-terminal branching. 

We first use the general structure of the ir, p-model to derive an expression for the probability 
of finding a tree 7 at a time t. We use this expression to prove that the mean m(n) and the 
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variance of (n) in the number of terminal segments at time t do not depend on ir s , provided p 7 is 
only dependent on the number of terminal segments and not on the topology of a tree. For this 
case we also derive differential equations for the probabilities p(n,t), i.e. the probability of finding 
an arbitrary tree with n terminal segments at time t. We apply these methods to the cBES-model 
and compare the mean /J.t(n) and the variance of (n) found to known results and approximations 
for the dBES-model. Then we develop efficient algorithms for the evaluation of ^-dependencies 
for the cBES-model, and use these to evaluate the topology distribution. 

B. Probability of finding a specific branching history 

Assuming that initially we start from the simplest tree, i.e. the tree with only one terminal 
segment fa, we want to calculate for a specific tree 7 with a specified branching history B the 
probability that it can be realized by the branching Poisson process and has not yet branched 
further. We write a branching sequence B as, 

B= (fa, 61 , ft, fc, ■■■,&), (8) 

where the bj indicate which terminal segment of the tree fa branched to obtain fa + \ as indicated 
in the following branch history diagram: 

fa $ fa h ■ ■ ■ V /3 n -l V fa- (9) 

Before we write down general expressions describing the probability of finding a branching sequence, 
we describe how we can find the probabilities for two short sequences. First a sequence B2 = 
(fa,b±, fa) leading from a tree with one terminal segment to a tree with two terminal segments 
without further branching before time £2 happens with probability 



p{B 2 ,t 2 ,h) 
ft 2 



/ 2 A felA e- Aft i,/9i*i J] e-'W* 2 -* 1 ^! 

hufii /V^e-V 42 -' 1 ^. (10) 
Jo 



The probability density used is the product of two parts: first the branching probability den- 
sity of the initial tree Af, li/ g 1 e _Ai 'i'' 3 i <1 , and second the simultaneous survival probability density 
(llsg/32 e~ As -' 9 2(' 2 ~' 1 )^ of the final tree's segments from the moment of the first branching to the 
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end time. Notice that the redistribution of probability between the terminal segments of the final 
tree does not influence the final probability for finding the tree. Next we determine the probability 
of finding the branch sequence £>2-s>3 = (#2, &2, /%) from a tree fa with two terminal segments to a 
tree fa with three terminal segments over the time period t\ to t 3 , 

P(#2->3,*3,*l) = 

/* 3 A 62A e- Ab 2./32(*2-*i) 

X e -^,/3 2 (*2-tl) 

= A b2/?2 [ t3 e- p ^-^ e - p ^ t3 - t2) dt 2 . 
Jti 

(11) 

In the expression above we again find the survival probability density for the tree fa from the last 
branching event to the end exp(— J2 S £/3 3 ^s,p 3 (t — £2)) and the terminal segment 62 branching prob- 
ability density Ab 2)( g 2 exp(— Ab 2)( g 2 (t2 — £].))■ I n addition we need to include the survival probability 
density exp(— J2s€/3 2 ,s^b2 A s ,/3 2 (*2 — h)) for the non-branching segments of the tree fa. Integrating 
this probability density over all intermediate branching times gives us the full transition probabil- 
ity. In a similar vein we can now express the probability B 3 = (fa, b±, fa, 6 2 , fa) that the branching 
sequence had already taken place at the time £3 provided we had the tree fa at time to = 0, by 
taking the integral over the probability density for branching of the initial segment multiplied by 
the probability that branching of terminal segment 62 will lead to the tree fa, 

p(B 3 ,t 3 ,t ) = 

= [ tS X buh e-P^- t ^p(B 2 ^ 3 ,t 3 ,t 1 )dt 1 
J to 

= r 3 tftiA 6li/3l e-^i-*o) 
Jt 

Jti 

(12) 

We can apply the same type of construction to larger tree histories as well. Using p(B,j,t n ,t n -j) 
for the probability that tree /3 n -j+i develops during the time interval (t n -j,t n ) into the tree fa, 
we obtain 



p(B,j,t n ,t n -j) = 
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\ f Pfl ... {tn—j + 1 tn-j) 

x j - l,*n,*n-j+l)dt n -j+l, 

(13) 

and by setting p(23, 1, t n , i n -i) to the survival probability of the final tree, 

p(B,l,t n ,t n . l ) = e~ p ^ {tn - t "- 1 \ (14) 

the recursion is correctly closed. From these recursion formulas we can see that the 7^, ft depen- 
dencies can easily be split off by gathering them in what we call in anticipation of their role in the 
cBES-model centrifugal order factors 0(B,j), 

0(B,j)= ( J] n bi A . (15) 
\i=n-j+l J 

We will make this split to facilitate our analysis of the 7r-dependence of tree probabilities and 
expectation values. Splitting off the centrifugal order factor also leads to the introduction of a 
new recursively defined object: the p-dependent integral factor I(B, t) only depending on p s . 's and 
related to p(B,j, t n ,t n -j) by the relation 

p(B,j,t n ,t n -j) = 

0(B,j)I(B,j,t n ,t n ^). (16) 

The integral factor I(B,n,t n ,to) can now be found through solving the following recursion rules 
derived from the recursion defined in equations [13] and [HI 



n-jj — 

Pn — j-\-l (^n — j' + l ~tn—j ) 



I(B, j, t n ,t r 

Pn-j+l j e 

x 1(23, j - l,t n ,t n -j +1 )dt n - j+ i, (17) 

where we introduced p n _ 1+ i to abbreviate p„ . The recursion is terminated by the following 
condition, 

I{B, 1, i n , i n _i) = P (B, 1, t n , i n _i). (18) 

The notation used so far contains the bare minimum for the manipulation of the integrals in the 
recursion relations. As we are, however, mainly interested in probabilities at time t, we will use 
a reduced notation and write p(B, t) = p(B, n, t n = t,to = 0) for the probability of finding a tree 
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with history B at time t. The integral I(B,n,t n ,to) is independent of the ir s and only depends 
on the In the cBES-model these p 7 's depend solely on the number of terminal segments n 7 . 
Instead of writing I(B, t) we will therefore use I(n, t) for I(B, n, t n = t,to = 0) when dealing with 
models for which the p^s depend exclusively on the number of terminal segments. The centrifugal 
order factors are time-independent by definition and therefore contain no reference to time. If we 
use the full branching history we will write 0(B) for 0(B, n) dropping the explicit reference to the 
number of included branching steps. 

C. Exclusively terminal segment number dependent whole tree branching rates 

In the cBES-model and other models for which the whole tree branching rates p's depend 
exclusively on the number of terminal segments, the probability p(n, t) that an arbitrary tree at 
time t has n terminal segments is by construction independent of it, because the total branch 
rate does not depend on tree topology. The sum of all probabilities contributing to p(n, t) should 
therefore also be 7r-independent: 

P(<M) = E f(7,*) 

7|n 7 =n 

= '(M) E E 0(B). (19) 

7|n 7 =n BdH-f 

Here we introduced p(j, t) as the total probability of finding the tree 7 at time t. We can express 
p(j, t) as the sum over the set containing all branching histories leading to the tree 7, 

p( 7 ,t)= E P(*M): (2°) 

and we used the factorization into a 7r-dependent and a 7r-independent part from equation [161 
Furthermore, we know that if the p's depend exclusively on the number of terminal segments, all 
branching histories leading to n terminals use the same sequence of p values and therefore they 
lead to the same value of I(B,t) which we therefore write as I(n,t). To check directly that p(n,t) 
is independent from tt we evaluate the double sum over centrifugal order factors 0(B). To achieve 
this we first change from summing over trees and their branch sequences to summing over terminal 
segments. Observe that every branching sequence of n — 1 branch events occurs once, so we can 
replace the double sum by a single sum over branch sequences containing n — 1 branch events after 
which a tree 73 with n terminal segments results, 

E E 0(B) = E 0(B) (21) 

7|n 7 =n BdH-f B\n{^rs)=n 
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This sum can be replaced by the sum over all branching sequences that are one branch event shorter 
but with every contribution multiplied by the sum over the branch rates of the terminal segments 
of the tree resulting from the shorter branch sequence. If we repeat this step until we are left with 
the root segment we see that what remains is a 7r-independent and equal to 1, 

E °( B ) 

B\n(-y B )=n 

(71-2 
n **i,Pi e f-.A.-i 
3=1 s£/3„_i 

E 0(B) 

B\n(-fi3)=n— 1 

= 1 (22) 
Putting this back into equation [19] we find 

p(n,t)=I(n,t). (23) 
and therefore combining this with equation [16] we can now express p(^y,t) simply as, 

p( 7> t) = I(n,t) 0{B). (24) 

This relation states that in the cBES-model and other models for which the p's depend exclusively 
on the number of terminal segments, the probability of finding a tree is a fixed part of the total 
probability of finding trees with the same number of terminal segments. In other words if n(j) = 
n{^') then the ratio p{"f,t)/p(-y',t) is time independent. 

A further direct consequence of the above is that expectation values of functions / which are 
exclusively dependent on the number of terminal segments of trees are also 7r-independent, 

7 

n=oo 

= £ f(n)I(n,t). (25) 

n=l 

This has the important consequence that moments n(n x ) of n are ^-independent and therefore 
contain no information about the redistribution of the total branching rate over the terminal 
segments. 

For numerical evaluation we will need the dynamics of the probability pin, t) of finding a tree 
in the subpopulation of trees with n terminal segments. Under the assumption used above the rate 
at which trees leave this subpopulation is independent of the particular tree and proportional to p n 
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and the size of the subpopulation p(n, t)N popu i a tion with N popu i a tion equal to the total population 
size. We find, therefore, that the rate at which trees leave the n and enter the n + 1 terminal 
segment subpopulation is equal to p n p{n : t)M popu iation- On the other hand the rate at which trees 
join the n and leave the n — 1 terminal segment subpopulation is given by p n -ip(n — l,t)J\f popu i at ion- 
Combining these two rates and dividing out the total population size M popu iation, we find that the 
rate of change of p(n, t) is given by: 



with (n > 1). To suit our situation these equations need to be combined with appropriate initial 
conditions. As before we assume that initially all trees have a single terminal segment p(l,0) = 1, 
p(n > 1, t) = and that the total number of trees is conserved, which we ensure by setting p = 
and p(0, t) = to prevent us from creating new trees. 



Here we will show that for the cBES-model we can obtain exact solutions for E = or E = 1 
irrespective of the value of S. In the appendix we show that the dBES-model can be solved for 
these values of E aswell, but only when 5 = 0. 

We use the rate equations for the p(n,i)'s to rewrite the time derivative of the mean and the 
variance, 



dp(n, t) 



= p n ~ip(n - 1, t) - p n p(n, t) 



(26) 



D. Exactly solvable cases of terminal segment number development 



dp(n) 



dt 




n=oo 



n{p n -ip{n - 1) - p n p(n)) 



n=l 

n=oo n=oo 



^2 (n + l)p n p(n) - Y n PnP{n) 



n=0 n=l 



n=oo 



PnPin) = p(pn) 



n=l 



(27) 



and similarly, 



dp,(n 2 ) 
dt 
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= n2 (pn-ip{n - 1) - p n p{n)) 

n=l 

n=oo n=oo 

= ^ (n + l) 2 pnP(n) - n 2 p n p{n) 



n=0 n=l 

n=oo 



= ]T (2n + l)p n p(n) 

n=l 

= 2p,(np n ) + 

(28) 

and making the proper subtractions we get a differential equation for the temporal development 
of the variance, 

da 2 (n) d 2 2 
-Ai = - (M n 2 )-^(n)) 

= 2fi(np n ) + - 2p(n)p{p n ). 



(29) 



For the cBES-model we obtain, 

' //,! " ; - hp 



dt 

= 2bp(n 2 - E ) + bp(n 1 - E ) 
-2bfi(n)p(n 1 ~ E ). 



(30) 



In general this is not a closed system of equations but for E = or E = 1 this system of equations 
actually becomes a closed system. For E = 1 we obtain, 

dp{n) 



dt 
da 2 in) 
dt 



= b, 

= b, (31) 



with simple linear solutions, 



For £ = 0we obtain, 



p(n,t) = p + bt, 
a 2 (n,t) = ol+bt. (32) 



dp(n) 

— = bp(n), 

= 2ba 2 (n) + bp(n), (33) 
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with exponential solutions, 

H(n,t) = /x e w , 

a 2 (n,t) = (a% + fx )e m - // e W . (34) 



E. Evaluating it- dependence 

In our treatment until now we have given an explicit label to each terminal segment and thus 
have treated our trees as labeled trees Harding lf| . For the description of the history of branching 
dendritic trees this is correct as the individual terminal segments can be followed in time. For our 
evaluation of the 7r-dependence we will, however, also use unlabeled trees. We will do this for two 
reasons. The first follows from the typical data at which this analysis is aimed, i.e. reconstructed 
dendritic trees, which have no natural labeling of the terminal segments associated with them. 
The second reason is the strong reduction in the number of trees we obtain by using the unlabeled 
trees in our algorithm achieving a strong reduction in computer memory usage. We will use the 
unlabeled trees as equivalence classes on our labeled trees, i.e. we will say that those labeled trees 
which are related to each other by rotating the subtrees around branching points have the same 
topology, which we then describe by an unlabeled tree. To make the distinction clear we will use a 
tilde to denote such equivalence classes, i.e. 7 denotes an equivalence class or unlabeled tree and 
7 denotes a labeled tree. 

The summed centrifugal order factor 0(7) for a labeled tree 7 is given by, 

0( 7 ) = £ O(B). (35) 

We would like to calculate this quantity without explicitly constructing all histories. That is, we 
shall express it in the summed centrifugal factors of those trees that are directly preceding 7 in the 
histories leading up to 7. We introduce a new notation to show the feasibility of this idea. We 
write 7^(7) for the set of direct predecessors of a tree 7, and for ifbjfij we write (which for 

the cBES-model is a valid notation because under the terminal growth hypothesis two subsequent 
trees implicitly fix the terminal segment that branched). We use B~ l to denote the branching 
sequence with the last branch event removed. This allows us to rewrite the previous expression as 



Oil) = E *7lA^°(S 



-1\ 



BeH-, 



E E ^n-MB- 1 ) 



/3„-l= 7 ' 
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time(1/b) time(1/b) time(1/b) 




FIG. 1: Temporal development in the cBES-Model. Top row: temporal development for the prob- 
abilities p(n,t) for n — 1, 6, 11, 101 and three different values of E: E = (left), E = i (middle) and 
E = 1 (right) in all cases lower n- values lead to earlier p(n, i)-peaks, while late peaking traces might peak 
outside the window shown. Middle row: expectation value and variance for the number of terminal seg- 
ments for three different values of E (matching those in the top row) calculated using the first 1000 p(n, t)'s 
while keeping Phigh < 10~ 6 . Bottom row left: comparison expectation value (solid grey lines) with mean 
field prediction (dashed black lines) for five different E values: E = 0, 0.25, 0.5, 0.75, 1; at i^-values 0, 1 the 
mean field solutions coincide with the exact solution. Bottom row right: comparison of mean field solu- 
tion and numerical results; after an initial growth of the relative error for intermediate values of E , i.e. 
E = 0.25, 0.5, 0.75, the relative error attenuates. The standard deviation and mean for E = and E = 1 
correspond within the numerical error with the exact solutions. 
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= E E o^ 1 ) 

7 'G-p( 7 ) B£H 7 

/9n- 1=7' 

= Y, ^h'°^')- (36) 

7'eP(7) 

To obtain the last line above we made use of the fact that every history leading up to a predecessor 
tree 7' appears only once in H^. 

To find the predecessors of the tree 7 we express it in terms of its two subtrees 7; and j r as 
7 = (7/,7 r ). Now the set of direct predecessors P(7) of the tree 7 can be found from the direct 
predecessors of its subtrees as 

P(7)=P((-yi,7r)) = 



U (P>7r) U (J ( 7 /,p) 

VPeP( 7i ) / \pG^(7r) 



(37) 



provided we can find 7' = (ji,p) and 7' = (/?, 7 r ). In other words, we now need a method to find 
the information for a tree 7' from the knowledge we stored about its two subtrees. At this point the 
fact that we will store and calculate the information based on unlabeled trees becomes important. 
Before we continue with further details of the actual calculation, we therefore need to introduce 
some relevant properties of unlabeled trees. 

n 

We use an enumeration of unlabeled trees devised by Harding 16] and first used for dendrites by 
van Pelt and Verwer [391). T his enumeration corresponds with reverse lexicographical ordering (van 
Elburg and van Ooyen [34]) and is the natural order in the recursive tree construction process that 
we use in our algorithm. In this recursive construction process we construct trees with n terminals 
from trees with a lower number of terminals, by first combining trees with n — 1 terminals with 
trees with 1 terminal, then combining trees with n — 2 terminals with trees with 2-terminals and so 
on. We stop this process after combining trees with (n + l)/2 terminals with those with (n — l)/2 
terminals for an odd n and combining all unique pairs of trees with n/2 terminals for an even n. 
Where the trees with the largest number of terminal segments are traversed in the outer loop and 
the trees with the smaller number of terminal segments in the inner loop. If in these loops we 
traverse the trees in their own order of construction then the order of construction of the larger 
trees is entirely fixed. There is one important detail, when the subtrees have an equal number 
of terminals, we avoid constructing the same topology twice by combining trees in the outer loop 
with trees with a higher index in the inner loop. An unlabeled tree 7 is now fully characterized 
by its number of terminal segments n 7 and its position in the construction order of trees with the 
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same number of terminal segments. 

With the ordering above the index of an unlabeled tree 7 with n = can be related to the 
indices of its subtrees 7z,7 r with 7/ < 7 r by the following relation, 

n 

7 = 7r+ N 3 N n-j 
j=n t +l 

f (2N ni - „ 

+°n u n r < ^ 7Z 

+(l-5 nu n r )(li-l)N nr , (38) 

where we identified the unlabeled tree with its index, i.e. 7 denotes both an unlabeled tree and 
its index. We use N x to denote the number of unlabeled trees with x terminal segments, 5 to 
denote the Kronecker delta function and n\ to denote n 7i and n r to denote n^ r . We are now ready 
to address the calculation of the centrifugal order factor, because equation [38] allows us to find 
for a tree 7 all the topologies of its direct predecessors from the topologies of its subtrees direct 
predecessors. 

We start by adapting the construction process of the unlabeled trees to set up the necessary 
book-keeping. We shall use equation [36] to find the centrifugal order factor for a tree 7. For the 
calculation of nyy = Ttbrf = 2~ Sut > JCy we need to know the centrifugal order of all the terminal 
segments in 7' for the calculation of the normalization factor Cy = J2bey 2~ Sub , and we need to 
have the centrifugal order of the branching terminal segment Ub available separately. Now if we 
have the tree 7 = (7/,7 r ) then the list of centrifugal orders of the set of its terminal segments is 
related those of its subtrees as 

[vb]be~/ = Wb + l]f>e 7! + tyb + l]be>, (39) 

where we used square brackets to indicate the list of centrifugal orders, addition within square 
brackets to indicate a change in values in the list, and addition outside the brackets to indicate the 
joining of the two lists, while keeping all elements. In words, the centrifugal order of the terminal 
segments of a tree corresponds to the centrifugal order of its subtrees increased by one, which is 
the extra distance added by the new root element. For the calculation of the normalization factor 
the order of these centrifugal orders is immaterial and the same holds for the calculation of factors 
■Kjtji. Furthermore, for each labeled tree corresponding to the same topology this list will contain 
the same numbers with the same multiplicities. Therefore, we can store all necessary information 
on the basis of topology, i.e. we can do the bookkeeping on the basis of the unlabeled trees as 

Wbh ■= H + % + H + %• ( 4 o) 
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Similarly we store for each unlabeled tree its direct predecessors and the centrifugal order of 
the terminal segment whose branching connects the direct predecessor to the current tree in a list 
as 

i(vb lW ),l%>€V(~/)- ( 41 ) 

Admittedly, our notation becomes somewhat unwieldy here, but we need just one more line of it. 
Provided lists like these are stored for the constituting subtrees 7/, j r , we can calculate the list for 
a tree 7 using the following, 

[( t/ 6 T |y))7 / ] 7 'eP(7) = 

[(%i V + 1 )(7',7r-))] 7 'eP( 7i ) 

(42) 

Again, we build these lists only once for every topology, i.e. we do the book-keeping on the basis of 
unlabeled trees. We can now give the necessary list assignment in an informal shorthand notation 
which hides some details which should be clear from the previous equation: 

[(f&,7% := 

[(^ + l,(7',7r))] 7 , 

+[fo + l,(7l,7'))]v (43) 

When we execute this assignment we use the index relation given by equation [38] to replace the 
specification of the direct predecessor trees in the indices of its subtrees, e.g. (72,7'), by the index 
of the topology of the direct predecessor itself. 

In our current implementation we calculate the normalization constants and TTyy using the list 
mentioned before in which we stored information on centrifugal orders and direct predecessors. 
When given a new value of S, our algorithm starts with calculating summed centrifugal order 
factors for small tree topologies and then works towards topologies of larger and larger trees. To 
achieve this the topologies are visited in the original construction order of the unlabeled trees. 
When we arrive at a specific unlabeled tree we have the summed centrifugal order factors of all 
its predecessors at our disposal. We also know the predecessors and the centrifugal orders of the 
terminal segments that branched, and we can therefore calculate the centrifugal order factor for 
the current topology using equation [36l 

The summed centrifugal order factors we calculated represent the probability of finding a single 
representative labeled tree corresponding to a topology. The probability of finding a topology is 
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therefore given by the product of the summed centrifugal order factor with the number of labeled 
trees corresponding to this topology. 



Branch sequences (each a 

/">/"\T»-J-1 nildflATl AT I / I 1 1 1 1 

CUIlLlIlLldLlUIl (Jl 1 r £\ _L , 1 J ( 


mail) 


hah) 


3(2(1, l),l)-> 4(3(2(1, l),l),ir 
3(2(1,1), 4(3(1, 2(1, 1)),1) 
3(1,2(1,1)) -> 4(1, 3(2(1, 1),1)) 
3(1, 2(1,1)) -+4(1, 3(1, 2(1,1))) 
3(1, 2(1,1)) -> 4(2(1,1), 2(1, 
3(2(1,1), 1)^4(2(1,1), 2(1, 


m = 4 


h = 1 


m = 1 


h=2 



TABLE I: Possible histories for trees with 4 terminal segments and the resulting multiplicities and number 
of histories. In the first column branch sequences are indicated and the resulting trees marked with 7 denote 
the form which we use to denote unlabeled trees. The second column indicates for the unlabeled tree in 
the preceding column the multiplicity m a (7), i.e. the number of different labeled trees corresponding to 
the indicated unlabeled tree. The third column indicates the number of histories h a (^y) for a labeled tree 
equivalent to the unlabeled tree. The number of histories for an unlabeled tree is equal to the product of 
m a (j) and ^(7). 



Let us now revisit the topic of unlabeled and labeled trees briefly and illustrate some of the 
terminology with an appropriate notational device. We can describe the branching structure of 
a tree by specifying at each segment starting from the root the number of terminal segments it 
carries. The tree structure is captured by putting the segments in the subtrees in brackets after 
the segment from which they bifurcate. For example, a tree with a single terminal segment is 
simply denoted as 1. The unique tree with 2 terminal segments is denoted as 2(1,1). The labeled 
trees with 3 terminal segments are 3(2(1, 1),1) and 3(1,2(1,1))), but these correspond to a single 
unlabeled tree, because the only difference is in the order of the two subtrees. 

Table U shows all possible histories leading to trees with 4 terminal segments. This table also 
shows that the unlabeled tree 4(3(2(1, 1),1),1) has 3 alternative forms, i.e. there are 4 labeled 
trees corresponding to this tree. For a general tree 7 the multiplicity 711,(7), i- e - t ne number of 
equivalent trees (including the tree itself), is given by 772(7) = 2 U ( " where 2/(7) denotes the number 



of unbalanced nodes (van Pelt and Verwer [39j), meaning nodes at which the two subtrees have 
different topologies. Alternatively, the multiplicity of a topology 7 with subtopologies ji and tv 
springing from the root segment is given by, 



771(7) = m(7/)7n(7 r )2^ ^i*). 



(44) 
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The tree 4(2(1,1), 2(1,1)) has a multiplicity of 1 but it is the first tree which has more than one 
possible history. In general the number of alternative histories can be calculated from the rela- 



tion Harding 



This relation is a combinatorial consequence of the fact that when we combine histories of two 
subtrees there are (n 7; + n 7r — 2)! possible orders in which we can select branching events from 
the two histories. However, we are not free to choose the order in which we pick these events from 
the two histories because in a branching history the order of events does matter. To correct for 
this we need to divide this factor by the factors giving the number of forbidden permutations on 
the subtree histories, i.e. we need to divide it by ((n 7! — l)!(n 7r — 1)!). As the branching of the 
root segment is always the first event to take place and it introduces no combinatorial factor, the 
number of histories of a tree can simply be calculated from the product of the number of histories 
of the constituting subtrees and this combinatorial factor. 

These relations allow us to calculate the probability p(j) of finding a topology 7 as 

p(7) = m( 7 )0( 7 ). (46) 

We furthermore obtain two usefull constraints which can act as checks on our code, 

E P(7) = l, (47) 

j\n{ ; y)=n 

and 

E ™(7)M7) = (« - 1)!, (48) 

7|ra(7)=ra 

where ^(7) is the number of histories of a single labeled tree corresponding to the topology 7. 

A further acceleration of our calculations can be achieved using the following recursive relations 
for normalization constants and 7r 7 | 7 /'s: 

C 7 = 

6G7 



e s (c~ n + C 7r ), 



(49) 
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and 



e in 



= —cZ 



r 



P yrn. 

5, 



S(y b . ,+i) 



The impact of using these recursive relations on memory usage is probably limited because the 
storing of centrifugal order is replaced by the storing the same number of factors 7r 7 |y . As at present 
the limiting factor in our analysis of ^-dependence lies in memory usage, we did not implement 
this improvement in our current code. 

The code for this paper is available from ModelDB at 
http://senselab.med.yale.edu/modeldb via accession number 129071. 



III. NUMERICAL RESULTS 



Although the main results of this paper are mathematical and algorithmical in nature, we 
developed and implemented all the tools needed for the numerical evaluation. This section describes 
these. Because there is a natural split between E-dependency and S-dependency we will discuss 
these results separately. The B-dependency is not discussed as it can be adjusted by simply 
making another choice of time units. In all results discussed below we have therefore used 6 = 1 
for simplicity. 



A. E-dependencies 

For low values of n we can numerically solve the system of differential equations that we derived 
for the p(n, t)'s: 

dp(n,t) . , 

— — — = p n -ip(n- l,t) -p n p(n,t), (51) 

with < n < M. The highest value of n for which we can compute p(n, t) is determined by the 
available memory. This restricts the possibilities to calculate the mean and variance of the number 
of terminal segments when we try to evaluate longer development times. In particular we are 
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FIG. 2: S-dependence cBES-Model. Distribution of probability over the trees for different values of the 
S parameter and for different numbers of terminal segments. Top row: graphs for eleven terminal segments 
(N = II). Middle Row: graphs for seventeen terminal segments (N = 17). Bottom row: graphs for 
twenty three terminal segments (N = 23). Left column: distribution of probability over different summed 
electrotonic path lengths. Middle column: distribution of probability over the logarithm of multiplicity or 
equivalently number of unbalanced branch points. Right column: distribution of probabilities over the tree 
asymmetry index. In all graphs S varies from —1 to 1 with steps of 0.25. The S = — 1 distribution is not 
offset, but all other distributions have offsets which are increased in steps of 0.08. The visible difference in 
smoothness between the graphs in the different rows is partly due to the number of trees for the different N 
values, which are 207, 24631 and 3626149 for N = 11, N = 17 and N = 23, respectively. In all graphs the 
S = 0.5 distribution is indicated with filled symbols to act as an extra landmark. 
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limited by the probability of finding n- values outside the range of values used during the numerical 
integration of the differential equation, i.e. ouphi g h{t) = J2n>M v(. n ->t)-> with M the highest ra-value 
used. For the cBES-model the time at which Phi g h(t) reaches a predetermined threshold is mainly 
determined by E. This is because given a range of ra-values the stiffness of the set of equations 
varies strongly with E: if we choose M = 1000 the rates p n span three orders of magnitude at 
E = while they are all equal for E = 1. The stiffness of the E = system severely limits 
the possibility of lowering Phi g h(t) by enlarging the system of differential equations. We used the 
odel5s method of Matlab (The Mathworks Inc., Natick, Massachusetts, USA) to integrate the stiff 
differential equations for the first 1000 p(n,t)'s. We gathered the probability Phigh(t) i n a separate 
variable, and kept this leak to higher n- values below 10 -6 . 

Let us first look at the development of the probabilities p(n, t) of finding trees with different 
numbers of terminal segments. For different values of E these probabilities are shown in the top 
row of figure [TJ For E = shown on the left we see a very fast spread of probability from small trees 
to trees with more terminal segments. In fact we see that for the values shown the probability of 
finding a tree with e.g. 6 terminal segments peaks before it crosses (if it ever crosses) the probability 
of finding a tree with 1 segment. For all values of E shown the peak probabilities are reached later 
for higher numbers of terminal segments. If we compare the different graphs we also see that as 
expected the E = case spreads fastest, followed by E = 0.5, while the slowest spreading case 
shown is E = 1. If we look at peak heights we see that in the E = 1 case a large proportion of 
the probability is concentrated in the peaking component and we can expect a small spread in the 
number of terminal segments. In the E = these peaks are much lower and we see that many 
components are present with similar probabilities, which is indicative of a large spread. 

In section III Dl we found exact results with which we can verify our numerical results for E = 
and E = 1. For the numerically found mean values the exact expressions differ from the numerical 
results by maximally 0.01% for the E = case and less than 10~ 13 % for the E = 1 case. Similarly 
for the standard deviations the maximal differences are smaller than 0.2% and 10~ 4 %, respectively. 
In principle these errors can be further lowered by including extra p(n, t) values in the recursion, 
but in the E = case we expect that an exponentially growing number of p(n, t)'s needs to be 
calculated for every time interval At added. We expect that the numerical errors in the mean for 
the other values of E between E = and E = 1 are between the values just mentioned, so smaller 
than 0.01% and larger than 10 -13 %. The second row of figure [I] shows the actual development of 
the mean number of terminal segments /i(n) and its standard deviation <r(n), and indeed we see 
the effects expected on the basis of the first row of figures. Furthermore we see that except for the 
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case of E = the standard deviation grows slower than the mean. 

The errors in the mean number of terminal segments found above are so small that the numerical 
solution can be used to evaluate the quality of the mean field approximation proposed by van 
Pelt van Pelt and Uylings [37] in the context of the dBES-model. They suggested replacing fj,(n 1 ~ E ) 
with n(n) 1 ~ E in equation [30] leading to 

= bm l ~ E , (52) 

where we used a tilde to indicate the use of the mean field approximation. This equation is solved 
for E ^ 0, fjLQ > and t > by 

fi(n) = (bt + ji E ) l / E . (53) 

The bottom left graph of figure HJ shows a comparison between the mean field solutions and the 
numerical solution, showing a visible discrepancy between the two, larger than the expected error 
found in the numerical solution. Despite this visible difference, there is still reasonable agreement 
between the mean field solution and the numerical solution. This agreement is further illustrated 
in the bottom right graph which shows the ratio between the mean field solution and the numerical 
solution. From this graph we see that for the values shown the relative error first increases and 
then peaks at values below 10 % before it decreases. The graphs in the middle row of figure [T] 
show that for E = 0.5 the standard deviation grows slower than the mean. This observation might 
explain why we see the mean field solution improve over time. 



B. S-dependencies 

To start this section we restate one of our mathematical results. For a branching mechanism 
for which the overall branching rate of a tree depends solely on the number of terminal segments, 
two trees 7 and 7' with an equal number of terminal segments 71(7) = n{^') are found in a time 
independent ratio, i.e for this situation p('y,t)/p('y' ,t) is time independent. Therefore if the cBES- 
model is a correct description the S-dependence can be evaluated at any point in time at which 
sufficient number of trees of a given number of terminal segments are available. 

In the left column of figure [2] the distributions are plotted against summed electrotonic path 
length. To calculate an electrotonic path length we assumed a length of l s = At for each segment 
including intermediate segments, set the electrotonic length constant for the terminal segments all 
to the same value At and then calculated the electrotonic length constants for the other segments 
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related to the terminal segment electrotonic length by applying a branch power of R = 3/2, 

X(s) = \ t \ln{ S yl R , (54) 

where n(s) denotes the number of terminal segments in the subtrees connected to the segment s. 
This equation incorporates both Rail's power law (Rail and the dependence of the electrotonic 
length constant on dendritic diameter. Rail's power law relates the diameter of a dendritic segment 
d(s) to the diameter of its daughter segments d(s\) and d(s2) through d(s) R = d(si) R + d(s2) R - 



The electrotonic length constant is given by A = \J dr m / 4r where d is the dendritic radius , r m the 
specific membrane resistance and r a the axial resistance. This gives every segment a dimensionless 
electrotonic length of 

As = T7~T = n{8)-V™. (55) 
A(s) 

These dimensionless lengths are added together to a quantity P St for all the segments in the path 
from and including the root segment to and including the terminal segment St- Finally all the P St 's 
are added together to form the summed electrotonic path length 

SEP^) = Ps t - (56) 

■st€7 

In general a symmetric tree obeying Rail's branch power law will have a larger summed electrotonic 
path length because the average dendritic segment is much thinner leading to a large electrical 
separation between terminal segments and soma. With fixed segment lengths, terminal segment 
diameter d and membrane resistance r m the actual values are not material to the comparison 

pare 



between trees. We use the summed electrotonic path length to enable us to consistently com; 
different morphologies in an electrophysiologically relevant way (van Elburg and van Ooyen 
without claiming to use realistic electrotonic lengths. The summed electrotonic path length values 
were binned in 35 equal size bins between the maximum and the minimum values found for a given 
number of terminal segments. 

In the middle column of figure [2] the distributions are plotted against the number of unbalanced 
branch points or equivalently the log of the multiplicity of the labeled tree corresponding to a single 
topology. Here more symmetrical trees will lead to a low number of unbalanced branch points and 
hence lower values correspond to higher symmetry. For the multiplicity we see all possible values 
and no binning was applied. 

In the right column of figure [2] the distributions are plotted against tree asymmetry index (van 



Pelt and Schierwagen 



351 ]. van Pelt et al. 40(]). The tree asymmetry index is given by taking the 
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mean of the partition asymmetries over all intermediate segments Sj in the tree 7: 



A( 7 ) 



1 



(57) 



n(7) - 1 



Where the partition asymmetry A Si for an intermediate segment Sj is given by comparing the 
number of terminal segments in the two subtrees trees 7/ and 7 r through the expression, 



The values of the tree asymmetry index range from zero for perfectly symmetrical trees to close to 
one for the most asymmetrical trees. The tree asymmetry values were binned in 20 equal size bins 
between the maximum and the minimum values found for a given number of terminal segments. 
For the multiplicity we show all possible values and no binning was necessary. 

Figure [2] shows that when we increase the number of terminal segments from N = 11 in the 
top row to N = 23 in the bottom row the distributions become smoother. This is due to the 
strong increase in the number of trees with increasing terminal segment number. At S = all trees 
with the same number of terminal segments have the same probability and these curves therefore 
indicate how the population of labeled trees is distributed over the underlying parameter if they 
appear with equal probability. For N = 11 it is clearly visible that the population is not smoothly 
distributed over the summed electrotonic path length and the tree asymmetry index. For the tree 
asymmetry index even at N = 17 the distribution is not smooth, despite the presence of 24631 
different topologies. 

When we move from S = — 1 to 5 = 1, i.e. from a situation in which more distal terminal 
segments show higher branch rates to a situation where more proximal terminal segments show 
higher branch rates, we see higher summed electrotonic path lengths, lower numbers of unbalanced 
branch points and a lower tree asymmetry index. Although the three indicators we used in this 
section differ in their detailed ordering of trees, these indicators do agree on the ordering of pairs 
of trees which are widely separated and they all implement some notion of relative symmetry in 
the topology of tree. (The code provided with this paper produces plots which allow readers to 
compare these measures.) Thus, despite the local differences in the ordering of the trees, the results 
in figure [2] show very clearly that increasing S decreases tree asymmetry as measured by the three 
different indicators introduced above. 




(58) 
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IV. DISCUSSION 

To describe neurite branching processes in continuous time seems more natural because in the 
underlying biology there is no intrinsic time step. However, in evaluating the biological plausibility 
of discrete time versus continuous time formulations we need to be more careful. In the continuous 
time formulation information about a branching event is regarded as spreading instantaneously 
over the whole tree. Whether this is a valid approximation is decided by the differences between 
the relevant time scales. In principle some time passes before the increased demand for resources 
at recently branched terminal segments influences the availability of resources at other terminal 
segments. If the delays in resource availability correspond roughly with the time step used in the 
discrete time model, then the discrete time model can capture some of the delay effects. This 
requires however that the discrete time step is treated as a separate variable of the model. If 
the delays in resource availability are short compared to the average interbranching time interval 
then the continuous time formulation seems better fit for modeling branching. We can have a 
closer look at this problem by comparing B with the time scale of potential rate-limiting processes. 



Reported values for B in dendritic branching (e.g. van Pelt and Uylings [36], van Pelt et al. 411 ]) 
are in the order of 1 to 10 branching events per day. If, for example, we assume that diffusional 
processes, which are generally considered to be slow, are rate limiting, then we can give some 
order of magnitude estimates about the time involved in the spreading of resource limitations. We 
first follow Hely et al. fis| . they considered MAP2 to be rate limiting and used a slow diffusion 
constant of 1 um 2 s~ 1 for MAP2. When we use this diffusion constant in a simple dimensional 

n 

analysis argument (Hentschel and Fine |19|]) typical timescales needed to spread a change in MAP2 
availability over a tree of size 100 \xm found are on the order of t = (1 um 2 s~ 1 )(100 /im) 2 = 10 4 s ~ 
3 hours. In this scenario the discrete time model seems more appropriate. Even more so if we 
consider applying the model to larger apical dendrites of cortical pyramidal neurons for which the 
typical timescale by the same argument would be of the order of a full day. If on the other hand 



we follow Hentschel and Fine [19| and assume a diffusion constant close to that of calcium to be 
rate limiting and a similar dendritic size, then typical timescales are on the order of t = 10s. And 
we would be lead to conclude that the continuous time model would be a better approximation. 
As many of the underlying factors are at present unknown ultimately new experiments and more 
detailed biophysical models should decide in which situations the different formulations are better. 
It is, however, important to note that to our knowledge nobody has sofar singled out a conceptual 
role for the size of the time step in the dBES-model. 
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Furthermore, the value and temporal development of the basal branching rate b are important 
for comparison of our results with experimental data. Kliemann 25| showed how time dependence 



of B can be modeled by modeling the branching process as a Galton- Watson process in a varying 
environment, but did not include dependence on terminal segment number or centrifugal order. 
Fortunately, in the cBES-model such a time-varying b(t) has no influence on the structure of the 
model dynamics. In fact it is possible to solve the model assuming 6=1 and then calculate the 
effective time T as the integrated basal branching rate T = Jq b(t)dt and use T instead of t in the 
final expressions to capture time-varying basal branching rate. The reason that this is possible 
is that in all our expressions t appears in combination with b, and bt = f b(t)dt for constant b. 
This shows that we can take our integrals with respect to b(t) instead of t. For the dBES-model 



36| and they found that to 



temporal development of b(t) has been studied by van Pelt and Uylings 
fit the experimental data a 'rapidly and monotonically decreasing function of time' is necessary. We 
expect that this experimental result carries over to the cBES-model with perhaps minor parameter 
changes. 

Terminal segment branching without pruning is the main case we analyzed here. A large part 
of our treatment, however, also applies to intermediate segment branching. Counting the number 
of potential histories of a certain topology will be complicated by intermediate terminal branching, 
but will still be tractable and one would be able to evaluate ^-dependencies or 7r-dependencies. 
The dependence on terminal segment number or p would be unchanged and can be evaluated in 

3e possible to analyze models 
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the same way as we showed for the cBES-model. It would thus 
like the centrifugal order dependent model by van Pelt and Verwer 
branching. Including the pruning of terminal segments would pose serious problems to our analysis. 
Evaluating the centrifugal order dependencies in the probability of the realization of a specific tree, 
for example, would need a summation over an infinite number of histories and progress in this area 
would require control over these infinite sums. In the presence of pruning the evaluation of functions 
depending solely on the number of terminal segments under the assumption that p{^) = p(rij) is 
still possible using the techniques presented here. But extra assumptions will have to be made 
about the creation and destruction of unbranched root segments. Including pruning does open the 
possibility of a stable distribution without a rapidly and monotonically decreasing basal branch 
rate. An initial investigation seems to indicate that in special cases equilibrium solutions for the 
distribution of the number of terminal segments can be found using detailed balance. We are, 
however, not aware of all the limitations of approaches using detailed balance and a thorough 
knowledge of these limitations is necessary to carry out such an analysis to the full. 
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As indicated above experiments and more detailed biophysical models are needed for a compar- 
ison of the continuous and discrete time formulations in specific contexts. However, before such 
comparisons can be made statistical tools need to be developed and/or implemented. Although the 
development of such statistical tools is outside the scope of this paper we think that the results pre- 
sented here are an important prerequisite. Furthermore, we think that the work presented here can 
contribute to the further development of reliable neural network simulators based on stochastically 
generated single cell morphologies. 

Appendix A: Details of the 7r, p-model 
1. Exact I(B,t) and p{B,t) 

We start this subsection with the derivation of explicit expressions for I(B,t) and p(B,t), for 
two special cases, the first case being that of a constant p n = b and the second case assuming that 
all p n for a particular branching sequence B are different, i.e. pi / pj if i 7^ j . These explicit 
expressions are given here because they played an important role in the inception phase of the 
work presented in this paper and because they have a usefull structure which can facilitate future 
work. The general case, which we will not discuss here because the cBES-model always gives rise 
to one of the special cases discussed here, will be a complicated mixture of these two results. The 
numerical methods we used for solving these systems of differential equations, although subject to 
other limitations, are not sensitive to the essential difference underlying the exact analysis of the 
two special cases. 

We know from section III CI that I(n,t) = p(n,t). For constant p n = b (corresponding to the 
E = 1 case in the cBES-model), we furthermore know that the distribution of branch events 
in time is the same as encountered in an ordinary Poisson process, i.e. the number of branch 
events during a time period is Poisson-distributed and therefore the number of terminal segments 
is Poisson-distributed: 

I(n,t) =p(n,i) 

where p n = b for all n. For the case where pi 7^ pj if i 7^ j, a condition which applies to the 
cBES-model provided E 7^ 1, we can derive explicit expressions for I(B,t) as well. All integrals we 
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encounter are of the following type: 

rU, 

in - 



e 



Pn — j-\-l {tn—j-\-l tn—j ) 



x e- p * ( - tn - tn -3+ l) dt n ^ j+1 (A2) 

with i > n — j + 1 . Evaluating this integral yields 

S(B,i,n-j + l) 

x ^e _Pn_J ' +1 ^ n ~* n_J ^ — e~ Pi ^ tn ~ tn ~ j ^ (A3) 

which is a sum over exponentials which are of the same type as the right factor in our integral, 
except with the integration variable i ra _j+i replaced with t n -j and further contains a factor 

S(B,i,j) = —?—. (A4) 
Pi ~ Pj 

The next integration, if needed, is completely analogous to the one shown here but with j replaced 
by j + 1. From this result we can distill the full result by making the following observations: both 
limits of integration yield a factor S(i,n—j + l); the upper limit of integration comes with a change 
in exponential, i.e. pi gets replaced by p n -j+i', the lower limit comes with an extra factor —1 and 
as mentioned before t n -j+i gets replaced with t n -j for both limits. If we keep contributions from 
upper and lower limits to I(n, t) separated, the n—1 integrals in I(n, t) will lead to 2 n ~ 1 exponential 
terms which are fully determined if we specify at which integrals we took the upper limit. If we 
use the indices i of the integration variables dti to denote at which integrations we took the upper 
limit, then the tuple u = (n, ...,0), containing in descending order all the upper limits used to 
arrive at the term and for technical reasons the opening and closing values u\ = n, U[ as t = 0, can 
be used to express this term in S(B,i,j) as follows 

n n (as) 

i<l(u) jWi>j>u i+1 

where S is defined as before except for S(B,i,0), which equals 

S(B,i,0) = e -«(*»-*°). (A6) 

If we insert factors pi from the recursion relation and sum over the set U(n — 1) which contains all 
possible sequences of upper limits over n — 1 integrals, then we obtain the full integral, 

mt) = (ii Pl ) ]r 

\i=l / ueW(n-l) 

x n (- 1 ) n (-s(B, Ui ,j)). 

i<l(u) j\ui>j>u i+1 

(A7) 
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2. I(B,t) for general 7r, p-models 

The systems of differential equations for p(n, t) used to find I(B, t) for those cases where p 7 
solely depends on the number of terminal segments n can be used for general w, />models as well. 
Remember that the I(B, t) where explicitly introduced as 7r-independent mathematical objects, 
but they have the additional property that they only depend on p values actually in the branching 
sequence B. This seems a rather trivial observation but it has the important consequence that if 
we can calculate I(B, t) for a w, />model A and we know that another tt, p-model B has the same 
p- values in the branching sequence B, then we can calculate I(B, t) for model B by using model 
A. In particular, we can define a model A by setting p^ = p^ , i.e. we can pick A such that its 
p-values depend solely on the number of terminal segments. Then, using model A, we have direct 
correspondence between I(n,t) and p(n,t), or equivalently between /(£>, t) and p(n,t). Therefore, 
we can calculate the p(n,t) for a suitable model A and use these as the values for I(B,t). The 
advantages of this approach to the calculation of I(B, t) and with it p(^y, t) are: we avoid the 
complicated calculation of I(n,t) from the equations I All and IA7I derived in appendix IA 11 we keep 
numerical precision because the explicit equations IA1I and IA7I rely on an increasing number of 
cancellations between small terms with increasing n, and the relation we found between I(n, t) 
and p(n, t) is a more generic one not limited to the special cases discussed in the first part of this 
section. 

Appendix B: Limitations to the analysis of the discrete time BES-model 

For comparison we also present here a novel analysis of the original discrete time BES-model 
using a transfer operator formulation. This will help pinpoint in this model the obstacles towards 
a more extensive exact treatment, which motivated us to study the continuous time model. 

1. The case of E = 0, S = 

We start with the simple case E = 0,S = 0,B=p for two reasons. Firstly we need to make 
clear the correspondences between model and notation. Secondly the textbook result (e.g. Dehling 
and Kalma [lo|). 



r \2 2 2,2 
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(Bl) 

on Galton- Watson processes is an important check on the more general result we will obtain later. 
This textbook result relates average and variance at time t to those at t — 1 and those after 
the first time step, where it was assumed that at t = there is only one terminal segment. A 
better interpretation without reference t = 1 is to interpret fi\ as the average offspring from a 
terminal segment during one time step and o\ as the variance therein. This interpretation allows 
one to specify different initial tree topology probability distributions. To stress this preferred 
interpretation we will write ji s and a s instead of fii and o\. 

In a Galton- Watson process the individuals in a population at time step t are the offspring of the 
individuals at time step t — 1. Furthermore, an individual's offspring is determined by a probability 
law which is independent of its history and the number of individuals in the population. Although 
when Watson conceived this problem, the main question was whether the population would die 
out, in our case we model dendritic branching without pruning and as a result the population of 
terminal segments will never die out. 

For the case at hand the branching probability is independent of the tree topology and time, 
and we simply have 

Pt(s,7)=P- (B2) 

The full probability law changes as folows: 

Pd = 0, 
Pc = 1 - P, 
Pb = p. 

(B3) 

The subscripts d, c, b stand for dying (pruning) , continuation and branching, respectively. We wish 
to point out that we allow an individual to be among its own offspring. This important deviation 

n n 

in comparison with static stochastic models (Devaud et al. [111 ]. Kliemann [25]) is necessary to com- 
plete the mapping of our problem to a Galton- Watson process and is crucial to our interpretation 
of the generation as the time step. 

The probability generating function g(x) for branching of one terminal segment during a time 
step becomes 

g{x) = Pd + Pc% + PbX 2 = (1 — p)x + px 2 , (B4) 
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From this we obtain 

Ms = sfO-) = 1 +P, 

O* = «/'(l) + </(l)-</(l) 2 

= 2p + (1 +p) - (1 + pf = p(l -p). 

(B5) 

Assuming fio = 1, cxq = we get the following for the average temporal development: 

t H = (l+p) t . (B6) 

2 

If we look at the quantity Vt = + - (^ s _i) A*t we can deduce from equation IB1I that it has a very 
simple time evolution, 

v t = (J%v t -i- (B7) 
With the equation above we can also obtain the variance at time t: 

2 °i 

a t = Vt 1 TT^* 

a ' 2 -0£-a4) 



^ s (/i s - 1) v 

= (l-p)((l+p)2*-l-(l+p)*-l). 

(B8) 

Our analysis sofar cannot be applied to the more general problem because it is explicitly assumed 
that the branching probabilities at every terminal are the same for all terminals during all time 
steps. In the next subsection we wish to focus on the case where we have a free choice of E and S. 

2. General case with explicit time dependence 

We will relate average and variation to functions of the distribution in the previous time step; 
for the case of S = E = we will find back equation IB 1[ The key step is to describe the probability 
of finding a tree 7 at time t in terms of transition probabilities from other trees 7', 

p t (7) = E T *(7>7')f*-i(7')- (B9) 
7' 

For our purposes it will be necessary to lump together some of the transition probabilities. To this 
end we also introduce the following notation: 

n 7 =n 

TtM) = ]T r t ( 7 ,Y). 

(BIO) 
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Let us return to the definition of the average and the variance of the number of terminal 
segments: 

7 

at (n) = J2(n<Y-fH(n)) 2 Pt(j). (Bll) 

7 

Using lB9l we can cast fit in the form, 

M«) =EE^KT')K-i(7'). (B12) 

n -y' 

If we change the summation order we get 

Vt{n) = ^fi t ,~ i '{n)pt-i{l'), (B13) 

7' 

where 

ft , 7 ( n ) = ^nT ( (n, 7 ') (B14) 

n 

is the average of n at time t provided the tree 7' is the only tree present at time t — 1. We can now 
easily calculate /it~(n) using the stochastic independence of branching events at different segments 
of the tree: 

Mt, 7 ( n ) = Y.VtA n ) 

= n 7 + £n~ B+1 . (B15) 

Similarly, using the stochastic independence of the variances of branching events at different seg- 
ments of the tree we can calculate fit^(n 2 ): 

s£7 se-y 

= X!(Pt(s,7) ~Pt(s,l)) 
+ [E(l+ft(«,7)) 

\sG7 

= n 2 i + Bn- E+1 + 2Bn- E + 2 



+B 2 n- 2E+2 1 



7 



C*7*2 



(B16) 
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where C 7 *2 is given by multiplying all the centrifugal orders in the definition of C 7 by a factor two: 

c^2 = Y. 2 ~ 2Sus - ( B17 ) 

Inserting these results into equation IB 131 for both f/,t(n), we find that the growth of /it equals a 
moment of the distribution: 

tH(n) - iH-i(n) = Bfi t ^(n- E+1 ). (B18) 

We see that there are two special values for which this relation closes on itself and no other moments 
of the distribution are needed. For £ = 0we get exponential growth: 

^(n) = (l + S)*/io(n), (B19) 

and for £ = lwe get linear growth: 

Ht(n) = iJto(n)+tB. (B20) 

From these equations we can immediately deduce that the S'-dependence does not influence 
the average growth during one time step, but depending on the values of E the change in the 
distribution might influence the growth in a next time step. It is also clear from these formulas 
that the variance and hence the distribution of probabilities after one time step do depend on S. 

For the variance we obtain 

a ti n ) = $30*t,7( n2 )-A*?(«))Pt-i(7) 

7 

-fH-l{n)lM-l{ri~ 
+B 2 of_ 1 (n-^ 1 ) 



R 2„ (-2E+2C~f*2 



(B21) 



This equation should yield the right Galton- Watson behavior if we choose S = E = 0. For S = 
we have 

lf = ^ (B22) 
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and we obtain 

(B23) 

which indeed corresponds with equation IB 11 because = /jl s = 1 + B and <rf = a 2 = (B — B 2 ). 
For S = 0, E = 1 we get 

<? 2 t(n) = ^W+A-itl) 

+2B(nt-i{n) - iH-\(n)nt-x{\)) 

+Bm- X {i) 

R 2„ I C-y*2 \ 
-B 

= ^iW + B-BVi^ 1 ). 

(B24) 

From this we see that at these parameter values the growth of the variance for a population of 
large trees will be nearly linear and for a population of small trees the growth of the variance will 
be less than linear. For general S we have 

< ^ < 1 (B25) 

and therefore for E = 1 the variance will have an upper bound cr 2 (n) = cr 2 (n) + Bt and a lower 
bound a?(n) = crg(n) + (B - B 2 )t. 
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